home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / cagd_lib / sbsp_int.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-14  |  11.1 KB  |  286 lines
  1. /******************************************************************************
  2. * SBsp-Int.c - Bspline surface interpolation/approximation schemes.          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Feb. 94.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11. #include "extra_fn.h"
  12.  
  13. /*****************************************************************************
  14. * DESCRIPTION:                                                               M
  15. * Given a set of points, PtList, computes a Bspline surface of order UOrder  M
  16. * by VOrder that interpolates or least square approximates the given set of  M
  17. * points.                                                                 M
  18. *   PtList is a NULL terminated array of linked lists of CagdPtStruct        M
  19. * structs.                                                                  M
  20. *   All linked lists in PtList must have the same length.             M
  21. *   U direction of surface is associated with array, V with the linked       M
  22. * lists.                                                            M
  23. *   The size of the control mesh of the resulting Bspline surface defaults   M
  24. * to the number of points in PtList (if SrfUSize = SrfVSize = 0).         M
  25. *   However, either numbers can smaller to yield a least square              M
  26. * approximation of the gievn data set.                                      M
  27. *   The created curve can be parametrized as specified by ParamType.         M
  28. *                                                                            *
  29. * PARAMETERS:                                                                M
  30. *   PtList:      A NULL terminating array of linked list of points.          M
  31. *   UOrder:      Of the to be created surface.                               M
  32. *   VOrder:      Of the to be created surface.                               M
  33. *   SrfUSize:    U size of the to be created surface. Must be at least as    M
  34. *                large as the array PtList.                         M
  35. *   SrfVSize:    V size of the to be created surface. Must be at least as    M
  36. *                large as the length of each list in PtList.                 M
  37. *   ParamType:    Type of parametrization.                                   M
  38. *                                                                            *
  39. * RETURN VALUE:                                                              M
  40. *   CagdSrfStruct *:   Constructed interpolating/approximating surface.      M
  41. *                                                                            *
  42. * KEYWORDS:                                                                  M
  43. *   BspSrfInterpPts, interpolation, least square approximation               M
  44. *****************************************************************************/
  45. CagdSrfStruct *BspSrfInterpPts(CagdPtStruct **PtList,
  46.                    int UOrder,
  47.                    int VOrder,
  48.                    int SrfUSize,
  49.                    int SrfVSize,
  50.                    CagdParametrizationType ParamType)
  51. {
  52.     int i, NumLinkedLists, NumPtsInList;
  53.     CagdPtStruct *Pt, **PtPtr;
  54.     CagdRType *UKV, *VKV, *PtUKnots, *PtVKnots, *AveKV, *RU, *RV, *R2;
  55.     CagdSrfStruct *Srf;
  56.     CagdCtlPtStruct
  57.     *CtlPt = NULL,
  58.     *CtlPtList = NULL;
  59.  
  60.     for (NumLinkedLists = 0, PtPtr = PtList;
  61.      *PtPtr != NULL;
  62.      NumLinkedLists++, PtPtr++);
  63.     for (NumPtsInList = 0, Pt = *PtList;
  64.      Pt != NULL;
  65.      NumPtsInList++, Pt = Pt -> Pnext);
  66.  
  67.     if (SrfUSize == 0)
  68.     SrfUSize = NumLinkedLists;
  69.     if (SrfVSize == 0)
  70.     SrfVSize = NumPtsInList;
  71.     if (NumLinkedLists < 3 ||
  72.     NumLinkedLists < UOrder ||
  73.     NumLinkedLists < SrfUSize ||
  74.     SrfUSize < UOrder ||
  75.     NumPtsInList < 3 ||
  76.     NumPtsInList < VOrder ||
  77.     NumPtsInList < SrfVSize ||
  78.     SrfVSize < VOrder)
  79.     return NULL;
  80.  
  81.     RU = PtUKnots =
  82.     (CagdRType *) IritMalloc(sizeof(CagdRType) * NumLinkedLists);
  83.     RV = PtVKnots =
  84.     (CagdRType *) IritMalloc(sizeof(CagdRType) * NumPtsInList);
  85.  
  86.     /* Compute parameter values at which surface will interpolate PtList. */
  87.     switch (ParamType) {
  88.     case CAGD_CHORD_LEN_PARAM:
  89.         case CAGD_CENTRIPETAL_PARAM:
  90.         /* No real support for chord length although we might be able to */
  91.         /* do something useful better that uniform parametrization.      */
  92.     case CAGD_UNIFORM_PARAM:
  93.     default:
  94.         for (i = 0; i < NumLinkedLists; i++)
  95.         *RU++ = ((RealType) i) / (NumLinkedLists - 1);
  96.         for (i = 0; i < NumPtsInList; i++)
  97.         *RV++ = ((RealType) i) / (NumPtsInList - 1);
  98.         break;
  99.     }
  100.  
  101.     /* Construct the two knot vectors of the Bspline surface. */
  102.     UKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * (SrfUSize + UOrder));
  103.     VKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * (SrfVSize + VOrder));
  104.  
  105.     AveKV = BspKnotAverage(PtUKnots, NumLinkedLists,
  106.                NumLinkedLists + UOrder - SrfUSize - 1);
  107.     BspKnotAffineTrans(AveKV, SrfUSize - UOrder + 2, PtUKnots[0] - AveKV[0],
  108.                (PtUKnots[NumLinkedLists - 1] - PtUKnots[0]) /
  109.                (AveKV[SrfUSize - UOrder + 1] - AveKV[0]));
  110.     for (i = 0, RU = UKV, R2 = AveKV; i < UOrder; i++)
  111.     *RU++ = *R2;
  112.     for (i = 0; i < SrfUSize - UOrder; i++)
  113.     *RU++ = *++R2;
  114.     for (i = 0, R2++; i < UOrder; i++)
  115.     *RU++ = *R2;
  116.     IritFree((VoidPtr) AveKV);
  117.  
  118.     AveKV = BspKnotAverage(PtVKnots, NumPtsInList,
  119.                NumPtsInList + VOrder - SrfVSize - 1);
  120.     BspKnotAffineTrans(AveKV, SrfVSize - VOrder + 2, PtVKnots[0] - AveKV[0],
  121.                (PtVKnots[NumPtsInList - 1] - PtVKnots[0]) /
  122.                (AveKV[SrfVSize - VOrder + 1] - AveKV[0]));
  123.     for (i = 0, RV = VKV, R2 = AveKV; i < VOrder; i++)
  124.     *RV++ = *R2;
  125.     for (i = 0; i < SrfVSize - VOrder; i++)
  126.     *RV++ = *++R2;
  127.     for (i = 0, R2++; i < VOrder; i++)
  128.     *RV++ = *R2;
  129.     IritFree((VoidPtr) AveKV);
  130.  
  131.     /* Convert to control points in a linear list */
  132.     for (PtPtr = PtList; *PtPtr != NULL; PtPtr++) {
  133.     for (Pt = *PtPtr, i = 0; Pt != NULL; Pt = Pt -> Pnext, i++) {
  134.         int j;
  135.  
  136.         if (CtlPtList == NULL)
  137.         CtlPtList = CtlPt = CagdCtlPtNew(CAGD_PT_E3_TYPE);
  138.         else {
  139.         CtlPt ->Pnext = CagdCtlPtNew(CAGD_PT_E3_TYPE);
  140.         CtlPt = CtlPt -> Pnext;
  141.         }
  142.         for (j = 0; j < 3; j++)
  143.         CtlPt -> Coords[j + 1] = Pt -> Pt[j];
  144.     }
  145.     if (i != NumPtsInList) {
  146.         CagdCtlPtFreeList(CtlPtList);
  147.  
  148.         IritFree((VoidPtr *) PtUKnots);
  149.         IritFree((VoidPtr *) PtVKnots);
  150.         IritFree((VoidPtr *) UKV);
  151.         IritFree((VoidPtr *) VKV);
  152.     
  153.         CAGD_FATAL_ERROR(CAGD_ERR_PT_OR_LEN_MISMATCH);
  154.         return NULL;
  155.     }
  156.     }
  157.  
  158.     Srf = BspSrfInterpolate(CtlPtList, NumLinkedLists, NumPtsInList,
  159.                 PtUKnots, PtVKnots, UKV, VKV,
  160.                 SrfUSize, SrfVSize, UOrder, VOrder);
  161.     CagdCtlPtFreeList(CtlPtList);
  162.  
  163.     IritFree((VoidPtr *) PtUKnots);
  164.     IritFree((VoidPtr *) PtVKnots);
  165.     IritFree((VoidPtr *) UKV);
  166.     IritFree((VoidPtr *) VKV);
  167.  
  168.     return Srf;
  169. }
  170.  
  171. /*****************************************************************************
  172. * DESCRIPTION:                                                               M
  173. * Given a set of points on a rectangular grid, PtList, parameter values the  M
  174. * surface should interpolate or approximate these grid points, U/VParams,    M
  175. * the expected two knot vectors of the surface, U/VKV, the expected lengths  M
  176. * U/VLength and orders U/VOrder of the Bspline surface, computes the Bspline M
  177. * surface's coefficients.                             M
  178. *   All points in PtList are assumed of the same type.                 M
  179. *                                                                            *
  180. * PARAMETERS:                                                                M
  181. *   PtList:     A long linked list (NumUPts * NumVPts) of points to          M
  182. *               interpolated or least square approximate.                    M
  183. *   NumUPts:    Number of points in PtList in the U direction.               M
  184. *   NumVPts:    Number of points in PtList in the V direction.               M
  185. *   UParams:    Parameter at which surface should interpolate or             M
  186. *               approxximate PtList in the U direction.                      M
  187. *   VParams:    Parameter at which surface should interpolate or             M
  188. *               approxximate PtList in the V direction.                      M
  189. *   UKV:        Requested knot vector form the surface in the U direction.   M
  190. *   VKV:        Requested knot vector form the surface in the V direction.   M
  191. *   ULength:    Requested length of control mesh of surface in U direction.  M
  192. *   VLength:    Requested length of control mesh of surface in V direction.  M
  193. *   UOrder:     Requested order of surface in U direction.                   M
  194. *   VOrder:     Requested order of surface in U direction.                   M
  195. *                                         *
  196. * RETURN VALUE:                                                              M
  197. *   CagdSrfStruct *:   Constructed interpolating/approximating surface.      M
  198. *                                                                            *
  199. * KEYWORDS:                                                                  M
  200. *   BspSrfInterpolate, interpolation, least square approximation             M
  201. *****************************************************************************/
  202. CagdSrfStruct *BspSrfInterpolate(CagdCtlPtStruct *PtList,
  203.                  int NumUPts,
  204.                  int NumVPts,
  205.                  CagdRType *UParams,
  206.                  CagdRType *VParams,
  207.                  CagdRType *UKV,
  208.                  CagdRType *VKV,
  209.                  int ULength,
  210.                  int VLength,
  211.                  int UOrder,
  212.                  int VOrder)
  213. {
  214.     CagdPointType
  215.     PtType = PtList -> PtType;
  216.     CagdBType
  217.     IsRational = CAGD_IS_RATIONAL_PT(PtType);
  218.     int i, j,
  219.     NumCoords = CAGD_NUM_OF_PT_COORD(PtType);
  220.     CagdCtlPtStruct *Pt;
  221.     CagdRType **SPoints;
  222.     CagdSrfStruct *Srf;
  223.     CagdCrvStruct **Crvs;
  224.  
  225.     /* Construct the Bspline surface and its two knot vectors. */
  226.     Srf = BspSrfNew(ULength, VLength, UOrder, VOrder, PtType);
  227.     SPoints = Srf -> Points;
  228.     CAGD_GEN_COPY(Srf -> UKnotVector, UKV,
  229.           (ULength + UOrder) * sizeof(CagdRType));
  230.     CAGD_GEN_COPY(Srf -> VKnotVector, VKV,
  231.           (VLength + VOrder) * sizeof(CagdRType));
  232.  
  233.     /* Interpolate the rows as set of curves. */
  234.     Crvs = (CagdCrvStruct **) IritMalloc(sizeof(CagdCrvStruct *) * NumVPts);
  235.     for (i = 0, Pt = PtList; i < NumVPts; i++) {
  236.     Crvs[i] = BspCrvInterpolate(Pt, NumUPts, UParams, UKV, ULength,
  237.                     UOrder, FALSE);
  238.  
  239.     for (j = 0; j < NumUPts; j++)             /* Advance to next row. */
  240.         Pt = Pt -> Pnext;
  241.     }
  242.  
  243.     /* Interpolate the columns from the curves of the rows. */
  244.     for (i = 0; i < ULength; i++) {
  245.     int k;
  246.     CagdCrvStruct *Crv;
  247.     CagdCtlPtStruct
  248.         *CtlPtList = NULL,
  249.         *CtlPt = NULL;
  250.     CagdRType **CPoints;
  251.  
  252.     for (j = 0; j < NumVPts; j++) {
  253.         CagdRType
  254.         **Points = Crvs[j] -> Points;
  255.  
  256.         if (CtlPtList == NULL)
  257.         CtlPtList = CtlPt = CagdCtlPtNew(CAGD_PT_E3_TYPE);
  258.         else {
  259.         CtlPt ->Pnext = CagdCtlPtNew(CAGD_PT_E3_TYPE);
  260.         CtlPt = CtlPt -> Pnext;
  261.         }
  262.  
  263.         for (k = !IsRational; k <= NumCoords; k++)
  264.         CtlPt -> Coords[k] = Points[k][i];
  265.     }
  266.  
  267.     /* Interpolate the column, copy to mesh, and free the curve. */
  268.     Crv = BspCrvInterpolate(CtlPtList, NumVPts, VParams, VKV, VLength,
  269.                 VOrder, FALSE);
  270.     CPoints = Crv -> Points;
  271.     CagdCtlPtFreeList(CtlPtList);
  272.  
  273.     for (j = 0; j < VLength; j++)
  274.         for (k = !IsRational; k <= NumCoords; k++)
  275.         SPoints[k][i + j * ULength] = CPoints[k][j];
  276.  
  277.     CagdCrvFree(Crv);
  278.     }
  279.  
  280.     for (i = 0; i < NumVPts; i++)
  281.     CagdCrvFree(Crvs[i]);
  282.     IritFree((VoidPtr) Crvs);
  283.  
  284.     return Srf;
  285. }
  286.